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ABSTRACT 

We apply a novel adaptive mesh refinement code, AMRVAC, to numerically investigate the 
various evolutionary phases in the interaction of a relativistic shell with its surrounding cold 
Interstellar Medium (ISM). We do this for both ID isotropic as well as full 2D jetlike fire- 
■^j- ' ball models. This is relevant for Gamma Ray Bursts, and we demonstrate that, thanks to the 

^T) [ AMR strategy, we resolve the internal structure of the shocked shell-ISM matter, which will 

leave its imprint on the GRB afterglow. We determine the deceleration from an initial Lorentz 
factor y = 100 up to the almost Newtonian y ~ 0(2) phase of the flow. We present axisym- 
metric 2D shell evolutions, with the 2D extent characterized by their initial opening angle. In 
such jetlike GRB models, we discuss the differences with the ID isotropic GRB equivalents. 
These are mainly due to thermally induced sideways expansions of both the shocked shell 
. and shocked ISM regions. We found that the propagating 2D ultrarelativistic shell does not 

accrete all the surrounding medium located within its initial opening angle. Part of this ISM 
matter gets pushed away laterally and forms a wide bow-shock configuration with swirling 
flow patterns trailing the thin shell. The resulting shell deceleration is quite different from 
that found in isotropic GRB models. As long as the lateral shell expansion is merely due to 
^ \ ballistic spreading of the shell, isotropic and 2D models agree perfectly. As thermally induced 

expansions eventually lead to significantly higher lateral speeds, the 2D shell interacts with 
comparably more ISM matter and decelerates earlier than its isotropic counterpart. 

\ Key words: Gamma Rays: Afterglow, Hydrodynamics, Theory - ISM: jets and outflows 

■ Galaxies: jets, ISM - methods: numerical, relativity, AMR 



1 INTRODUCTION 

Many high energy astrophysical phenomena involve relativistic 
flows and shocks. For example, relativistic flows are invoked to 
explain the observed properties of various com pact astrophysical 
objects ( lAronsl2004lFerrarill99ilCorbell2004h . Astrophysical rel- 
ativistic flows can reach a Lorentz factor of 2 - 10 in associatio n 
with jets from Seyfert and radio loud galaxies JPiner et alj|2003h . 
or even go up to Lorentz factors 10 2 - 10 3 for Gamma Ray Burst 
(GRB) scenarios dSari & Pirarj ll999: Soderberg & Ramirez-Ruizl 
2001; Mes zaroll200r3r In the last decade, continued development 
of numerical algorithms and the increase in computer power have 
allowed to significantly progress in high-resolution, hydrodynamic 
numerical simulations in both special and general relativity (see 
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Marti & Mflller 2003). The enormous time and length scale ranges 
associated with violent astrophysical phenomena in relativistic hy- 
drodynamics (RHD), make Adaptive Mesh Refinement (AMR) an 
important algorithmic ingredient for computationally affordable 
simulations. RHD numerical simulations, particularly when com- 
bined with AMR capabilities, can investigate many details of rela- 
tivistic flow regimes relevant for astrophysics. 

In this paper, we concentrate on relativistic dynamics in the 
fireball model for the afterglow phases of GRBs, in one and two 
dimensio nal simulations. S ince the follow-up detection of GRBs 
in X-ray {Co sta et al. 1997) and their afterglows a t longer wave- 
lengths dSahu et all 1 19971: IVan Paradiis et all 1 19971 : iGalama et ail 
1 19971: iFrail et alj|l997l : IPiro et al.lll998l), the cosmologica l origin 
of GRBs has been established (Metz geTl997l : IWiierslll997h . These 
detecti ons confirmed t he predictions from the fireball theoretical 
model jRhoadsl 1 19931 : iKatzl 1 1994 iMeszaros & Reesl 1 1991 Ivietril 
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1 1997b . In this model, a compact source releases a large amount 
of energy in a very short timescale, producing a fireball expand- 
ing with relativistic velocity. Its internal energy gets fully con- 
verted to kinetic energy, leading to a shell expanding with very 
high Lorentz factor. This cold shell continues to expand and in- 
teract with the circumburst medium, producing a relativistic shock- 
dominated evolution. As the shell sweeps up the matter, it begins 
to decelerate. Here, we investigate the details of such propagating 
relati vistic shells with the re lativistic hydrodynami cs code AMR- 
VAC ( iBergmans etai]l2004l) . The AMRVAC code dKeppens et ail 
120031) is here for the first time applied to the numerically challeng- 
ing regime of high Lorentz factor, and we therefore include a va- 
riety of test problems, demonstrating the robustness as well as the 
limitations of our computational strategy. 

largely 



Up till recently, ana lysis of GRB flows have 
been d one analytically ( Shem i & Piraj Il99d: Sari & Piran 
1 19951 ; iMeszaros & Reesl 1 1991 IChiang & Dermerl \\99% . 
combined with numerical approaches usually employing a 
Lagrangian code. These latter works mainly investigate spher- 
ically symmet ric GRB scena rios for obvious compu t ational 

convenience dPanaitescu et al.l 1 19971 : iKobavashi etail 1 19991 ; 
iKobavashi & Sari 2000). Recently, some analytical works 
started to inves ti gate t h e multidimensional jet structur e in GRBs 
dRhoadsl fl997l (l999i: Ipanaitescu & Meszaroj [l999l: ISari et alj 



19991: iKumar & Panaitescj 120001: IPanaitescu & Kumaj I200I : 



Cheng et all l200ll : lOren et alj |2004 iKumar & Granotl [20030 . 
and some numerical simulations emerged as well, but restricted 
to relatively low (o rder 25) Lorentz factor dGranot et all 1200 ll ; 
ICannizzo et alj|20o"4T) . Higher speeds were obtained in the numer- 
ical simulation of the propagation of an axisymmetric jet through 
a coll apsing rotating massive star, as investigated by lAlov et all 
(2000) to analyse the first phase of GRBs. In these simulations, the 
jet is further followed after breakout to a maximum Lorentz factor 
of 7max ~ 44, which is still relatively small to the values required 
for GRBs by the fireball model. Therefore, an important area of 
current investigations in GRB context is to model the dynamics of 
narrow jets of ultra-relativistically flowing ejecta. This is motivated 
by the need to reduce the total amount of energy released in GRBs, 
by assuming these jets to point towards the observer, as compared 
to fully isotropic equivalents. This n eed is particularly cl ear for 
the exemp lary cases of GRB 990123 dKulkami et aHll999l). GRB 
050820A dCenko et alj2006l) . and GRB 050904 dFrail et al.ll2006f) . 

The detection of polarization 1 Covino et al.ll 999: Wiie rs et"al] 
1 19991 : lGreineretal1l2003l : lLazzati et alj |2004|) gave further sup- 
port to the jetlike model. Evidence for narrow collimated outflows 
in GRBs is sustained also by the achromatic breaks in the after- 
glow light curves which was predicted analytically (Rhoads (1997, 
1999); Sa ri et al. (1999)) and then observed in a large number 
of GRBs dStanek et al.lll999l; ISari et alj 1 19991: Igerger et alj|2000l: 



IPanaitesc JT2005I ; 1 120061 : iBarthelmv et alj l2005h . In iBloom etall 
d2003l) , various GRBs were analysed and in 16 of them, the com- 
bination of these breaks in the spectrum and the jet-like model 
was used to deduce their effective energy, which was about E ~ 
10 51 ergs. The half opening angle of such jets in GRBs is inferred 
to be of order few degrees. As a result, the afterglow produc- 
ing shocke d region is collimated too, with a similar initial open- 
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ing angle (Frail et al 2001; Berger et al. 2003; Cenko et al. 2006; 



20051 ; iDar & De R uiula 20(3). In our 2D simulations, 



we will concentrate on the afterglow phases in the GRB evolution 
starting from collimated ejecta, and discuss those dynamical effects 
causing opening angle changes in detail. Direct comparison with 



the evolution of an equivalent ID spherical shell is enlightening in 
this respect. 

This paper is organised as follows. We start by reviewing the 
relativistic hydrodynamic equations. In Section [3] we include sev- 
eral tests to demonstrate the AMRVAC code potential for realistic 
RHD computations. In Section|4] we present our main astrophysi- 
cal application to GRB flows in ID and 2D models. 



2 RELATIVISTIC HYDRODYNAMIC EQUATIONS 

The special relativistic hydrodynamic evolution of a perfect fluid 
is governed by the conservation of the number of particles, and 
energy-momentum conservation. These two conservation laws can 
be written as 



(pu% = 0, (T" r ) M = 0. 



(1) 



where p, it = (y,yrf), and T 1 " = phuf u v + pg^' define, respec- 
tively, the proper density, the four-velocity and the stress-energy 
tensor of the perfect fluid. Their definition involves the Lorentz 
factor 7, the fluid pressure p, and the relativistic specific enthalpy 
h = 1 + e + pip where e is the specific internal energy. For the 
(inverse) metric gf v , we take the Minkowski metric. Units are taken 
where the light speed equals unity. 

These equations can be written in conservative form involv- 
ing the Cartesian coordinate axes and the time axis of a fixed 'lab' 
Lorentzian reference frame as 



dU dF> n 
— + > = . 

at axJ 



7=1 



(2) 



The conserved variables can be taken as 

U = = yp, 3 = y 2 p hv,T = y 2 p h - p - 7pj , 

and the fluxes are then given by 

F = [py y 2 p hvv + pi, y 2 p hv - yp i?J , 

where I is the 3x3 identity matrix. To close this system of equa- 
tions, we use the equation of state (EOS) for an ideal gas, which is 
the polytropic equation with the polytropic index T, 



(3) 



(4) 



p = (T-l)pe. 



(5) 



At each time step in the numerical integration, the primitive 
variables (p, v, p) involved in flux expressions should be derived 
from the conservative variables U resulting in a system of nonlinear 
equations. One can bring this system into a single equation for the 
pressure p, 

^ n Mn P + Tp(y(pf -i) 
t + D- y{p) D — j = ' W 

which, once solved for p yields i? = T+ S p+D ■ This nonlinear equa- 
tion {6]l is solved using a Newton-Raphson algorithm. 



3 TESTING AMRVAC 

In view of the challenges in the numerical investigation of rela- 
tivistic fluids, we include here several substantial test results for 
code validation. We performed a large series of tests, some of them 
shown in this section. An important subclass of test cases is formed 
by Riemann problems, whose numerical solution can be compared 
to analytical solutions. Other, 2D tests shown here have no known 
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Figure 1. One-dimensional relativistic shock problem in planar geometry 
at t = 0.36. The solid lines are the analytical solution. 



analytical solution. Therefore, we compare the results of our sim- 
ulations with similar results previously obtained by other codes as 
documented in the astrophysical literature. 

The Adaptive Mesh Refinement version of the Versatile Ad- 
vection Code (AMRVAC) is specifically designed for simulating 
dynamics governed by a system of (near-)conservation laws. Avail- 
able equations are the Euler and magnetohydrodynamic systems, 
in both classical and special relativistic versions. The discretiza- 
tion is finite volume based, and various shock-capturing algorithms 
can be used. The automated AMR strategies implemented vary 
from th e original patch-based to a nov el hybrid block-based ap- 
proach jVan der Hoist & Ke ppens 2006). These procedures gener- 
ate or destruct hierarchically nested grids with subsequently finer 
mesh spacing. The refinement criterio n used in AMRVAC i s based 
on a Richardson type error estimator ( Kep pens etal]|2003h . In all 
following tests, we use a refinement ratio of 2 between consecutive 
levels, unless stated otherwise. 



Table 1. L\ errors of the density for the ID Riemann problem 1 with uni- 
form grid shown at t = 0.36 



Number of grid points 


u 


200 


1.15 x 10-' 


400 


6.4 x 10~ 2 


800 


3.2 x 10~ 2 


1600 


1.9 x ur 2 


3200 


1.06 x 10~ 2 



nonvanishing tangential velocities, for two ideal gases with poly- 
tropic index T = 5/3. We separate two different constant states 
p L = 10 3 , p L = 1.0 (left) and p R = 10~ 2 , p R = 1.0 (right). For 
the transverse velocity we form n i ne combinations of th e pair v Vi l 
and Vvr. As in iPons et alj fcood) : iMignone et alj d2005t) . we take 
v y , L = (0.0, 0.9, 0.99)c in combination with v„, R = (0.0, 0.9, 0.99)c. 
The spatial separation between the two states is initially at x = 0.5. 
The results at t = 0.4 are shown in Fig. [2] where we also overplot 
the exact solution using the code in lMarti & Muller 1 20031). 

The relativistic effects in these tests are mainly thermodynam- 
ical in the first mild test, and are due to coupling between the ther- 
modynamics (through specific enthalpy) and kinetic properties (by 
the initial tangential velocities). For small tangential velocity cases, 
we use only a resolution of 200 cells on the base level and 4 levels. 
However, for a high tangential velocity case, we use high base res- 
olution 400 with 10 levels to resolve the contact discontinuity and 
the tail of the rarefaction wave. In fact, for a high tangential veloc- 
ity at left (in the high pressure state), the effective inertia of the left 
state increases. This makes the occurring shock move slower and 
decreases the distance between the tail of th e rarefaction wave and 
the co ntact discontinuity. As also found in Izhang & M acFadven 
(2006), it remains a numerical challenge to capture the contact dis- 
continuity properly, which we only managed here by allowing a 
very high effective resolution. 



3.1 One-dimensional test problems 

3.1.1 Riemann problems 

In ID Riemann problems, we follow the evolution of an initial dis- 
continuity between two constant thermodynamical states. In ID 
RHD, we then typically find the appearance of up to three non- 
linear waves. Generally, one finds a shock wave propagating into 
the lower density/pressure medium, a rarefaction wave propagat- 
ing at sound speed into the denser medium, and between these two 
states, there can be a contact discontinuity. In the tests that follow, 
calculations are done in Cartesian geometry on a spatial domain 
< x < 1. The exact solutions for Riemann problems in relativistic 
hydrodynamics are discussed for vanis hing tangential speed ( i.e. y 



Marti & Muller (1994) and 



Pons etal.1 (12000). 



or/and z components for velocities) in 
for arbitrary tangential flow velocity in J 

In a first, mild test, we assume an ideal gas with polytropic 
index 5/3 and initial constant states characterized by p L = 13.3, 
p L = 10.0 (left) and p R = 0.66 x 10~ 6 ,p R = 1.0 (right), separated 
at the location x = 0.5. The results at t = 0.36 are shown in Fig.Q] 
with a resolution of 100 cells on the base level and 4 levels, where 
we also overplot the exact solution. In the table 1, we present the 
L[ = X (A*;) \pj -p(xf)\ norm errors of the density p, where p(Xj) is 
th e exact solution. The accura c y of our result is comparable to that 
of iLucas-Serrano et alj j2004h ; lzhang & MacFadven] J2006h - 

In a second test, we look particularly into effects due to 



3.1.2 Shock Heating Test 

In another ID test case, a cold fluid hits a wall and a shock front 
propagates back into the fluid, compressing and heating it as the 
kinetic energy is converted into internal energy. Behind the shock, 
the fluid becomes at rest. This t est has an analytical solu tion in 
planar symmetry as considered by Blandfor d & McKeejdl97 6). and 
the jump conditions are 



Pi 



Pi(n - D(nr+i) 
nr + i 



Pr 



r- 1 



cr- i) 



TiVi 



(7) 



7i + 1 

These give the post shock pressure p 2 and density p 2 values in terms 
of the incoming density and Lorentz factor, together with the shock 
propagation velocity v s h- 

In our test we take the same init ial conditions as in the recent 
paper bv lZhang & MacFadvenl I2006), where a cold fluid p = 10~ 4 
with a density p = 1.0 has an impact velocity of vi = (l.O - 10" 10 ). 
This corresponds to a Lorentz factor y = 70710.675. The tem- 
perature after the shock becomes relativistic, and therefore we 
take the polytropic index T = 4/3. Hence the shock velocity is 
v s h = 0.33332862. The AMR simulation is done with 20 cells on 
the base level and 4 levels on the spatial range < x < 1. The 
result at t = 2, with the reflective wall at x = 1, is shown in Fig. [3] 
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The exact solution is overplotted as well. In this test, because of 
the constant state behind the shock, the maximum impact Lorentz 
factor that can be achieved numerically is limited only by the pre- 
cision of the Newton-Raphson subroutine. This test is important to 
demonstrate its accurate treatment, in view of the intended simula- 
tions aimed at afterglows in GRBs. Indeed, in the shell-frame, the 
circumburst medium hits the dense shell with a high Lorentz factor. 
In a process similar to what is found in the above test, the kinetic 
energy of the impacting medium is converted to thermal energy of 
the external medium. Viewed in the lab frame, the swept up circum- 
burst medium will have similarly high Lorentz factor and will form 
a hot shocked layer ahead of the contact interface. Note also that 
Fig.|3]indicates that our discretization and w all treatment does not 
suffer from the visible density errors seen in lZhang & MacFadvenl 
d2006h . 



3.2 Two-dimensional tests 

3.2.1 A relativistic 2D Riemann problem 

A two dimensional square region is divided into four equal areas 
with a constant state each. We fix the polytropic index T = 5/3 
and assume free outflow boundar y conditions. The relativis t ic ver - 
sion of this test was proposed by Del Zann a & Bucciantinil J2002I) 



and subsequently reproduc ed by iLucas-Serrano et al. J2004) : 
IZhang & M acFadven (2006J) a nd un der slightly improved initial 
conditions by iMignone et alj (2005J). We repeat this simulation 
with t he same initial configuration from lDel Zanna & Bucciantinil 
d2002h . namely 



(p,V x /c,v y /c,p) = (0.1,0.0,0.0,0.01), 

/ \NW 

(p, vjc, Vyjc, p) = (0.1, 0.99, 0.0, 1.0) , 



v/c 


P/8e9 


p/5 e5 







Figure 3. One-dimensional shock heating problem in planar geometry, 
where a cold fluid hits a wall located at x = 1 . The results presented corre- 
spond to t = 2. The computational grid consists of 20 zones with 4 levels of 
refinement. The solid lines are the analytical solution. 



[p,V X /C,Vy/C,pj 

{p,vjc,v y /c,pj 



(0.5,0.0,0.0,1.0), 
(0.1,0.0,0.99,1.0). 



(8) 



The simulation is done with 48 x 48 cells at the lowest grid 
level, and we allow for 4 levels. The result is shown in Fig. [4] Our 
result is in qualitative agreement with those results published, and 
shows the stationary contact discontinuities between SW-SE and 
SW-NW with a jump in the transverse velocity. These are some- 
what diffused by the employed Total Variation Diminishing Lax- 
Friedrichs (TVDLF) discretization (Toth & Odstrcil 2006). A sim- 
ple and easily affordable remedy for improvement is to activate 
many more grid levels. Shocks feature across the interfaces NW- 
NE and SE-NE, propagating diagonally to the NE region, and an 
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Figure 4. Density distribution for the two dimensional shock tube problem 
at t = 0.4. With a polytropic index T = 5/3, a base resolution of 48 X 48 and 
4 AMR levels. 



elongated diagonal shock structure forms as the NE sector recedes 
into the RHS top corner. In the SW corner, an oblique jet-like struc- 
ture forms with a bow shock. 



3.2.2 Relativistic jet in 2D cylindrical geometry 

Since it is relevant for our 2D GRB simulations, we also present 
a two-dimensional simulation of an axisymmetric relativistic jet 
propa gating in a unifor m medium. We simulate the C2 jet model 
from lMarti et alj (1997), but with an enlarged domain and at higher 
effective resolution. Our computational domain covers the region 
< r < 15 and < z < 50 jet radii. Initially, the relativis- 
tic jet beam occupies the region r < 1,2 < 1, with vj et = 0.99, 
Pj et = 0.01 and its classical Mach number M = 6. In this case, the 
jet is super-sonic but its temperature is still classical, so we can take 
the polytropic index T = 5/3. The density of the external medium 
is Pcxt = 1-0. We follow the evolution until / = 130, and this end 
result is shown in Fig. [5] We performed the simulation with a reso- 
lution at the lowest level of the grid set to 90 x 300, and allowed for 
a total of 5 levels of refinement eventually achieving an effective 
resolution of 1440 x 4800. 

In this simulation, the relativistic motion of the flow domi- 
nates, the thermal energy is weak compared to the kinetic energy. 
As a result the external medium influences only weakly the jet and 
the Lorentz factor y ~ 7 flow produces a cocoon structure from 
the tip of the jet. One also finds a weak tra nsverse expansion of the 
outflow in accord with what is reported bv lMarti et al.1 dl997l) . This 
transverse expans ion of the jet is induced b y the pressure build-up 
inside the cocoon Begelman & Cioffil ( 1989). In our simulation, the 
average transverse expansion obtained is occurring at an estimated 
speed of Vt = 0.11 c. Moreover, at the contact interface between 
shocked external medium and jet material, complex vortical struc- 
tures form. These originate from Kelvin-Helmoltz type instabili- 
ties, as a consequence of cold fast jet outflow meeting a more static 
medium. The average propagation speed of the jet head is found to 
be 0.414 c, which is in agreement with the one-dimen sional analyt- 
ical estimate of 0.42 c as given bv lMarti etaIUl997t) . 




min-0.021025, max-705.051453 



loglO(rho) 



Figure 5. Density distribution for the axisymmetric relativistic jet at t = 
130. At left, we show the lab frame density, at right, we show the proper 
density in a logarithmic scale. The computational base grid consists of 90 X 
300 zones with 5 levels of refinement and the domain size is 15 X 50. 



L °g(p) 




Figure 6. Density distribution, in logarithmic scale, for the forward facing 
step problem, at t = 4.26. 
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3.2.3 Wind tunnel with step 

We reproduce here a standard test in the hydro dynami c litera - 
ture, namely the forward f acing step test from lEmervl dl968h : 
IWoodward & Colellal dl984l). but adjusted to the relativistic hy- 
dro re gime as in lLucas-Serrano et al.l d2004l) : IZhang & MacFadven 
d2006h . A horizontal relativistic supersonic flow enters a tunnel 
with a flat forward facing step. The test was done with a resolution 
50 x 100 zones with 4 levels. The size of the tunnel is < x < 3 and 
< y < 1. The step is 0.2 in height and its position is at x = 0.6. It 
is treated as a reflecting boundary. The upper boundary and lower 
boundary for x < 0.6 are also both reflecting. However, the left 
boundary is fixed at the given inflow and the right one has free out- 
flow. Initially, the whole computational domain is filled with ideal 
gas with T = 7/5 with a density p = 1.4 moving at v v = 0.999, 
i.e., with a Lorentz factor y = 22.31. The Newtonian Mach num- 
ber is set to 3. The result of our simulation is shown in Fig. [6] at 
time t = 4.26. In this test, the relativistic flow collides with the 
step, as a result a reverse shock propagates back against the flow 
direction and this shock reflects from the upper boundary. A Mach 
stem forms and remains stationar y. The result of our simula tion is 
comparable to what is reported in lLucas-Serrano et al. 



4 GRBS AND MODELS FOR THEIR AFTERGLOW 
PHASE 

A popular model for GRB flows is known as the fireball model. In 
this model, GRBs are produced by a relativistic outflow following 
a violent event near a compact object. A large amount of energy 
is promptly released by the com pact source in a region with small 
baryon loading (for a review see lPiranl |2005)). Initially, most of the 
energy of the flow is in the form of internal (thermal) energy. The 
shell expands rapidly converting its internal energy to kinetic. After 
the acceleration phase is complete, the shell is cold and moves with 
relativistic speeds. 

This cold shell interacts with the circumburst medium, pro- 
ducing strong shocks. Our simulations will consider the dynamics 
from this phase onwards. As the shell sweeps up mass from the 
external medium, the kinetic energy in the relativistic shell is grad- 
ually transferred to kinetic and internal energy in the shocked am- 
bient medium. Moreover, the shell itself gets traversed by a reverse 
shock, which in turn converts the kinetic energy of the shell to in- 
ternal energy. 

The observed afterglow emission that follows the prompt GRB 
emission is believed to come from synchrotron (with possible in- 
verse Compton contribution) emitting electrons that are accelerated 
in the forward and reverse shocks dSari et ai][l9 98: Gala maet al.l 
1998). In the initial phases of the shell-ISM interaction, the elec- 
trons can be in the fast cooling regime (i.e. their cooling timescale 
is shorter than the expansion timescale) and, therefore, radiate effi- 
ciently most of the energy injected to them. Furthermore, if most of 
the energy dissipated in the shocks accelerates the electrons, then 
one has to consider radiative shocks. If either of the previous con- 
ditions does not hold, the radiative losses in the shocks are small. 
Here, we assume that the radiative losses are dynamically unim- 
portant, i,e., the shocks are adiabatic throughout these simulations. 
According to the magnetization of the shell, the int eraction shell- 
ISM a nd the spectrum could change as is shown in iMimica et al.l 
(2006). Here, we assume that the magnetic field is dynamically 
unimportant. 



4.1 ID isotropic shell evolution 

In this simulation, we consider an ISM with uniform number den- 
sity 71ism = 1 cirT 3 . Many GRB afterglows (more t han 25%) seem 
to be produced in such constant density medium {Chevalier & Lil 
l2000l : IPanaitescu & KumaJl2002l ; IChevalier et alj|2004h . This con- 
stant density medium can be the resultant of a Wolf-Rayet star 
progenitor, with its su rroundings shaped by a weak stellar wind 
dVan Marie etalj|2006l) . Initially we set a uniform relativistic shell 
at Rp = 10 16 cm fro m the central engine, since according to 
IWoods & Loebl d 19951) the interaction of the shell with the ISM be- 
comes appreciable at this distance. The shell has an initial Lorentz 
factor of 7 = 100 (a 7 > 10 is in accord with a shell which is 
optica lly thin to gamma-rays dWoods & Loebl 1 19951 : lsari & PirarJ 
1995)), and energy 

E = 10 54 ergs = 47r 7 2 F&p+m*. (9) 

where S stands for the lab-frame thickness of the shell set to 

5 x 10 12 cm is of the order of the expected value for a fireball 

6 ~ max(c At, Ro/y 2 ), where At is the duration of the GRB. The 
ISM and the shell are cold, and the initial pressure is set to pi SM = 
10~ 3 ;jism nip c 2 and p s h c ii = 10~ 3 rc s heii ni p c 2 respectively. Note that 
this implies a huge initial contrast in the density measured in the 
lab-frame between the shell and the ISM D s hdi/AsM ~ 10', an d 
this presents an extreme challenge from a computational point of 
view. Initially, the energy of the shell is then mainly kinetic. We use 
a constant polytropic index T = 4/3, as the interaction shell-ISM 
will be dominated by the forward shock, where the temperature of 
the shocked ISM becomes relativistic. 

In this simulation we use an effective resolution of 1536000 
cells corresponding to the highest grid level 10 allowed. We use 
the full AMR capabilities in this simulation, since we simulate on 
a domain of size [0.3, 300] x 10 16 cm, with 30000 grid points on the 
lowest level. At f = 0, the shell itself is then only resolved from 
grid level 6 onwards, when we use a refinement ratio of 2 between 
consecutive levels. The initial shell is resolved by about 25 cells in 
grid level 10 (later in the dynamical evolution this means that there 
are many more grid points throughout the widening structure). We 
use this very high effective resolution to avoid any numerical diffu- 
sion which may cause an artificial spreading of the shell. We ensure 
that throughout the entire simulation, grid level 10 is activated and 
concentrates fine grids on both the forward shock and reverse shock 
regions. Both are very important to determine the precise timing of 
the deceleration. 

In Fig. [7] we show snapshots taken at lab-frame time t = 

2.2 x 10 6 s corresponding to an early time in the entire simulation, 
and in Fig. [9] we show snapshots taken at time t 1.5 x 10 7 s cor- 
responding to a time when the shock is fully developed, we will 
concentrate our discussion mainly on this figure [9] These figures 
demonstrate that we resolve all four regions that characterise the 
interaction between an outward moving relativistic shell and the 
cold ISM. From right to left, we recognize (Fig. [9j (1) the ISM 
at rest, (2) the shocked ISM that has passed through the forward 
shock, with its Lorentz factor raised to y^ ~ 30. This swept-up 
ISM gets compressed at the front shock and its number density 
reaches n {2) ~ 75cnT 3 ~ ^gli dSari & Piradl 19951) . These two 
values correspond to the analytical estimate given by eq. (|7J for the 
front shock propagation. Region (3) represents part of the initial 
shell material which is shocked by the reverse shock. The reverse 
shock propagates back into the cold shell, reducing its Lorentz fac- 
tor and converting its kinetic to thermal energy. Transfer of energy 
from the initial cold shell thus occurs both at the forward and the re- 
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Figure 7. The five regions characterizing the interaction of the relativistic shell with the ISM at t = 2.2 X 10 6 s. Panel(a): Lorentz factor, (b): logarithm of the 
pressure, (c): logarithm of density. 
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Figure 8. Ratio of thermal energy to mass energy, when the shell reaches distance: R ^ 6.5 X 10 16 cm (; ^ 2.2 X l() 6 s) (left), and R =s 4.3 X 10 17 cm 
(f=* 1.5 X 10 7 s) (right). 



verse s hock. Regions (2) and (3) are separated by a contact discon- 
tinuity (M eszaros &Reej| 19921) . At this spherical contact surface, 
the longitudinal velocity (Lorentz factor in ID case) and pressure 
remain constant, but there is a jump in density. Furthermore, re- 
gion (4) is the unshocked cold material of the shell, moving with 
a Lorentz factor 7 (4l = 100. The weak thermal energy of the shell 
interior itself did induce a slight expansion in the thickness of this 
part of the shell. 

The reverse shock separating region (3) and (4) prop- 
agates into the cold shell with a Lorentz factor y RS = 
7(3)7(4) (l - v ( 3)V(4)/c 2 ) ~ 2.5. This reverse shock is Newtonian 
inefficient in raising the thermal energy content as is shown in 
Fig. [8] (left panel), where we draw the specific thermal energy 
in the shocked ISM and shell when the shell reaches a distance 
R ^ 6.5 x 10 16 cm. The reverse shock remains Newtonian until it 
reaches a distance from the GRB source of R ~ 3.8 x 10 17 cm. Then 
it becomes mildly relativistic until R ~ 4.3 x 10 17 cm where the 
reverse shock becomes very efficient to convert the kinetic energy 
to thermal energy (see Fig[8]at right). Beyond this latter distance, 
the density of the unshocked shell part p (4) has decreased in accord 
with the spherical expansion of the shell, to p (4 ) <K 7 2 4) Pism- As a 
result, the reverse shock becomes relativistic. This behavior is char- 
acteristic for an initial thin cold relativistic shell decelerating in a 
constant density external medium. In fact, until the outward prop- 
agating shell reaches R ~ 3.8 x 10 l7 cm, the shocked ISM matter 
is hot e ( 2> ~ 30.0p(2) (where 7(2) = 30 corresponds to the analyti- 
cal solution for the relativistic forward shock e (2 ) = (7(2) - 1) P(2))> 
while the shocked shell material which has e (3) = e (2) is cold, since 
e ( 3j ~ O.Olpo). When the density in the non-shocked shell (i.e. p (4 j) 



decreases enough due to spherical expansion, the Lorentz factor of 
the reverse shock increases and the last part of the shocked shell 
becomes hot e (3 ) ~ p (3) . 

There is another region (5) indicated in the figures behind the 
shell. The density and the pressure in the region (5) are very small 
with n(5)_ m i„ < 10~ 6 cnr 3 and P(5) -ra j„ < 10~ 6 m p c 2 . Therefore, the 
region (5) is in the numerical point of view a vacuum. This region 
(5) is not of strong interest for the physics of the afterglow, but 
it is computationally challenging to resolve the interface between 
the regions (4) and (5) where the ratio of the lab frame density 
between the two reaches Z)( 4 )/Z)(5) ~ 10 14 in the first phase of the 
propagation of the shell (R ~ 10 16 cm), while the expansion of the 
shell remains weak. By the time shown in Fig. [7] this contrast has 
dropped to a value of at most 10 12 . 

The near-total deceleration of the shell only takes place when 
the two shocks in the shell-ISM interaction manage to convert an 
important fraction of the kinetic energy of the shell to thermal 
energy (and the efficiency of this conversion depends on whether 
the reverse shock is relativistic or Newtonian, as discussed above), 
while the rest is transferred to the swept-up ISM in the form of ki- 
netic and thermal energy. In the first phase of the deceleration, the 
maximum Lorentz factor of the shell decreases gently from 100 at a 
distance of R ~ 2.5x 10 I7 cm to 80 at a distance of R ~ 4.3x 10 17 cm. 
However, only at the latter distance of 4.3 x 10 17 cm, a sudden de- 
crease of the maximum Lorentz factor of the entire configuration 
from 7 = 80 to 7 = 30 takes place. This fast drop of the maxi- 
mum Lorentz factor as seen in Fig. [TO] coincides with the moment 
at which the reverse shock reaches the back end of the cold shell, 
thereby converting its kinetic to thermal energy (Fig. [9}. In fact, 
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Figure 9. The five zones present when the relativistic shell interacts with the ISM at t =s 1.5 X 10 7 s. (a) Lorentz factor, (b) log of pressure, (c) log of density. 
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Figure 10. The variation of the maximal Lorentz factor in the propa- 
gating shell-ISM structure with time a) when the shell propagates from 
R = 10 16 cm to R = 300 X Rq b) when the shell decelerates following 
the Blandford-McKee profile. 



when plotting the maximal Lorentz factor as a function of distance, 
initially we always observe the Lorentz factor of the unshocked 
shell matter. As soon as the reverse shock has crossed the entire 
initial shell, we start to follow the evolution of the Lorentz factor 
of the shocked ISM (at the forward shock) where the maximum 
Lorentz factor is 30 at that particular moment. 



After this phase, the shell structure continues to decelerate by 
transferring its kinetic energy to shocked ISM matter at the forward 
shock. However, as seen in Fig.QJJ] it still takes a certain time be- 
fore the variation of maximum Lorentz factor now characterizing 
the shocked ISM matter follows the self-simi lar analytical solution 
for bla st-wave deceleration as put forward bv lBlandford & McKed 
(1976). From about a distance of 1.2 x 10 18 cm, our numerical solu- 
tion starts to follow the analytical solution precisely. In fact, af- 
ter the reverse shock traversed the entire initial shell, a forward 
traveling rarefaction wave propagates through the entire structure 
thereby slowing it down while transferring most of the energy 
to shocked ISM regions. This structure does not follow the self- 
similar prescription and causes the initial difference. In the end, the 
distance between the forward shock and the contact discontinuity 
increased sufficiently and the resulting radial thermodynamic pro- 
files in between become fully described by the Blandford-McKee 
analytical solution. The Lorentz factor predicted by the Blandford- 
McKee solution behind the forward shock (for an adiabatic shock) 
is Tbm = (E/pisM c2 R 3 ) [l2 K R^ 12 - The prediction of the Blandford- 
McKee solution for the Lorentz factor of the flow is also plotted in 
Fig. [j0]and the agreement with the results of the simulation at these 
later stages of the decelation is good. 

Eventually, we enter into the mildly relativistic regime for the 
blast wave evolution. The transition to the Sedov-Taylor phase oc- 
curs beyond the simulated distance R > 300x 10 16 cm, since we still 
have a Lorentz factor of about 3 at the end of the simulation. The 
Sedov-Taylor distance we find is close to the analytical estimate 
given by / ~ (3£/4wpi S M c 2 )" 3 ~ 5 x 10 18 cm. 



4.2 2D modeling of directed ejecta 

Precise analysis of the afterglow phases requires to evolve numer- 
ically confined ejecta in more than ID, propagating in a jet-like 
fashion into the ISM. We now present axisymmetric, 2D simula- 
tions of a relativistic cold shell propagating in uniform ISM with 
a number density w ISM = 1 cm -3 . In this work, we investigate the 
uniform model jet (Rhoads 1999). The shell density and energy is 
set constant throughout the shell, and we take it to correspond to 
an isotropic spherical shell containing an equivalent isotropic en- 
ergy E iso = 10 5l ergs and a Lorentz factor y = 100. To make the 
2D computation feasible, we now start the simulation with a shell 
thickness 5 = 10 14 cm at a distance R = 10 16 cm from the cen- 
tral engine. In the initial setup the shell occupies an annular region, 
with half opening angle of the shell equal to 6 = 1°. This angle is 
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rather small with respect to those typically deduced from model- 
ing of the optical light-curve breaks but still in agreement with th e 
most collimated GRB flows teloom et ai1l2003l : IPanaitescdl2005l) . 
With the choice of a rather narrow jet, we expect the 2D effects to 
appear earlier and to be more pronounced with respect to a spher- 
ical shell with the same isotropic equivalent energy. Note that as a 
result, this jet like outflow has a decreased effective energy in the 
shell E ja = (6 2 /2)E iso . 

The AMR run uses 4 grid levels, taking 400 x 6000 at level 1, 
but with refinement ratios of 2, 4 and 2 between consecutive levels 
eventually achieving an effective resolution of 6400 x 96000. The 
domain is [0,4] x [0.3,30] in R units, as we specifically intend 
to model in detail the most dramatic phase of deceleration prior 
to the Blanford-McKee evolution. Note that we initially have y > 
1 /&, which appears to be the case for a GRB jet. Beaming effects 
are invoked to explain the observed steepening in the decay light 
curve resultant from the transition y > 1/6 (indistinguishable from 
an isotropic explosion) to y < 1/ 6 (Panaitescu & Meszarosl l 19991 ; 
IPanaitescu & Kumaj|2002l ; |Panaitescij|2005 ). With these setup we 
can verify from our high resolution simulation whether up to times 
corresponding to the transition y ~ 1/6, only the isotropic e nergy 
£i so is relevant for the dynamics and the resulting emission ( IPiranl 
l200dlGranotll2005h . 

The initial velocity of the shell is purely radial. Note that, com- 
pared to the ID isotropic case presented in the previous section, the 
2D simulation starts with an initial condition containing less en- 
ergy. This is done for mere practical reasons: we wish to keep the 
computation feasible within two week's execution time on a single 
processor. Due to this lower energy content, the deceleration dis- 
tance will be smaller by about one order of magnitude since less 
swept-up ISM mass is sufficient to decelerate the shell. This re- 
duces the need for resolving many decades of propagation distance 
as measured in units of the initial shell thickness. In fact the ID 
equivalent isotropic case with the same energy shows a sudden de- 
crease of the maximum Lorentz factor of the decelerating configu- 
ration that corresponds to the reverse shock crossing of the shell at 
R ~ 9x 10' 7 cm. This happens well before the simulated 3x 10 17 cm. 

Fig. QT] shows at the (top), the sound speed contour, and the 
density distribution in a logarithmic scale in the (left), and the 
Lorentz factor in the (right). At the bottom, a zoom in the region 
around the shell is shown in the (left), the sound speed and the 
lateral velocity, in the (right) the Lorentz factor, and in the (cen- 
ter) a zoom only on the shell, the density and Lorentz factor. As 
in the ID case, at first the shell propagates with a constant maxi- 
mum Lorentz factor, and this is accompanied by a weak spread of 
the shell. Part of this radial shell widening in the bottom part of 
the shell is affected by the creation of a very low pressure and den- 
sity region below the shell (also occurring in the ID scenario). This 
near-vacuum state remains at the rear part as the shell moves away 
at the specified Lorentz factor. In this 2D simulation, the unshocked 
shell also spreads laterally with an initial transverse (horizontal) ve- 
locity, since the shell is launched with a pure radial velocity. The 
corresponding maximum initial lateral velocity of the unshocked 
shell is Vt = 0.0175 c. However, the shocked, swept up ISM matter 
spreads laterally much faster, due to its high thermal energy con- 
tent. Initially, that shocked ISM part spreads with a comoving ve- 
locity of Vt,co ~ 0.4c, which is less than the maximum sound speed 
allowed by the polytropic equation of state c/ V3. Due to this fast 
sideways expansion of shocked shell and ISM, the mass of the ISM 
hit by the shell grows faster than r 1 . Therefore, the deceleration of 
the shell starts earlier than in the isotropic case, see FigQ/2] This re- 
sult implies that the transition from the phase where E lso is relevant, 



to the phase where Ej St is relevant in the dynamics takes place when 
the shocked ISM and shell start to spread laterally much faster than 
what corresponds to pure radial (ballistic) flow. 

In the last part of the shell-ISM deceleration phase, when the 
reverse shock has crossed the entire initial shell material, the lat- 
eral velocity of shocked shell material reaches a comoving speed 
of Vt, co ~ 0.7c. This means that we do find that the lateral ve- 
locity can be bigger than the sound sp eed in the mediu m which 
is in accord with the analytical result of ISari et alj 11999 | ). This is 
at odd s with numerical findings as those found in lCannizzo et al.1 
( 2004), which employ a much reduced resolution as compared to 
our AMR results (at low resolution, we do obtain a reduced lat- 
eral spreading velocity). As a result of this fast lateral spread of 
the shocked material, distinct differences occur in the deceleration 
stages as compared to the isotropic case. This result is very impor- 
tant, as it shows that the lateral spreading of the shell is not related 
only to the Lorentz factor of the shell but to the type of the reverse 
shock. In our computation, in an early phase the reverse shock is 
Newtonian and the expansion of the shocked shell part is modest. 
However, in a later phase the reverse shock becomes relativistic and 
this leads to faster lateral spreads. However, as the forward shock 
is always relativistic, already in an early stage the shocked ISM 
spreads with high velocity. The overall spreading of the shocked 
ISM and shocked shell configuration can, thus, be quite complex 
and rather more evolved than the one that semi-analytical models 
dRhoadsll999MPanaitescu & Meszarosll999l ; ISari et aljl999k|plrar] 
l200d) predict 

Only that part of the ISM found within the solid angle of the 
expanding shell is swept up, and this opening angle changes due 
to the spreading effects just discussed. In our 2D simulation, we 
found in analogy with the ID (higher energy) case from above, that 
the reverse shock is initially Newtonian, so the thermal energy in 
the shocked part of the shell does not increase a lot and its lateral 
expansion remains weak for a while. Later on, its lateral expansion 
speed goes up to the 0.7c mentioned above, as the reverse shock 
becomes relativistic and the material through which the shocked 
shell expands laterally has already been brought to lower densities 
by the shocked ISM interaction. 

The variation of the maximum Lorentz factor is less sud- 
den than in the equivalent ID spherical explosion, as quantified in 
Fig. Q/2] The part of the shell most distant from the symmetry axis 
decelerates before the more internal part. The shell sweeps up more 
matter than in the corresponding isotropic case in the external parts 
due to the lateral spreading effects discussed. In fact, in this simu- 
lation we may draw the analogy between the shell interaction with 
the ISM and simulations of relativistic AGN jet propagation into 
an external medium. As in those cases, the energy is transferred 
to the ISM through a bow shock structure. However, our modeled 
ejected shell representative of a burst in GRBs is not continually 
supported by injection of energy at the bottom. As a result, in the 
(small opening angle) shell there is not really evidence of a clear 
jet beam as in an AGN jet. In this case the interaction shell-ISM 
is dominated primarily by forward, reverse shock pair, and contact 
discontinuity in between. The changing 2D structure of this shock 
leads to differences in the shocked ISM mass loaded on to the shell. 
As stated earlier, only a fraction of the ISM within the opening an- 
gle of the shell is swept up, and an important part of ISM matter 
gets pushed away laterally as the thin radially confined shell ad- 
vances. The resulting behavior is clearly influenced by these 2D 
effects, and is the reason why the maximum Lorentz factor of the 
configuration starts to decrease earlier than in the ID scenario. In 
an isotropic scenario all swept up mass of the ISM remains in front 



© 2007 RAS, MNRAS 000, ??-?? 



10 Z Meliani, R. Keppens, F. Casse and D. Giannios 





Figure 11. On the (top), the sound speed, logarithm of density (right), and Lorentz factor contours for the 2D simulation (left). On the (bottom), a zoom 
around on the relativistic shell, the sound speed and the lateral velocity (left), the Lorentz factor (right), and the zoom on the shell, the density and Lorentz 
factor (center) at t = 5 x 10 6 s. 
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Figure 12. The variation of the maximum Lorentz factor evolution in ID 
and 2D scenarios compared. 



of the initial shell where most of the energy of the shell is continu- 
ally transferred to shocked ISM. However, in the jet-like explosion, 
the shocked ISM and shell expand laterally leading to interaction 
with more ISM material, but the main part of this ISM material gets 
deflected about the shell. As seen in FigQT| we find rather complex 
flow patterns trailing the thin shell. Hence, the mass accreted on 
to the shell in fact decreases as compared to the equivalent local 
isotropic scenario. Therefore, less energy is transferred continually 
to the shocked ISM. 

One interesting characteristic of the maximum Lorentz factor 
of the beamed shell at radii R > 7 x 10 16 cm is the "bumps" that it 
shows as function of radius. These modulations of y max are a result 
of rapid internal motions of the decelerating configuration caused 
by the complex shell-ISM interaction. In view to the very rich and 
unexpected early afterglow phenomenology re vealed by the SWIFT 
satellite (see, for example. IZhang et alj ( l2006h ). it is interesting to 
study whether these proper motions can cause a significant modu- 
lation in the emitted radiation expected from these flows. 

In this simulation, we find no real indication of a strong 
change in lateral spreading of the shell whe n the Lorentz 
factor drops down to y c = 1/6 ~ 57 tehoadsl 1 19991 ; 
IPanaitescu & Meszaros 199§). In fact, an important change is pro- 
duced later, when the Lorentz factor of the shell becomes smaller 
than 30, while the lateral velocity of the shell reaches vt. co ~ 0.7c. 
As pointed out, this coincides with the time when the reverse shock 
became relativistic and almost crossed the initial shell entirely. Af- 
ter this phase of rapid lateral spread, we find that the spread out 
shell decelerates faster, as it accumulates more matter (Fig. II U . 
More simulations with different values for the shell opening an- 
gle and thickness are to investigate how these parameters affect the 
phase of the deceleration where the lateral spreading of the shell 
becomes important. 



5 CONCLUSIONS 

In this paper, we presented and applied the AMRVAC code in its 
extension to relativistic hydrodynamics. The adaptive mesh refine- 
ment is particularly useful for simulating highly relativistic flow dy- 
namics. We always used the robust TVDLF sheme, and this shock- 
capturing method together with high effective resolution delivers 
numerical results that can rival or even improve other high order 



methods. As is well-known, difficulties in special relativistic hy- 
drodynamic simulations result from the non linear coupling be- 
tween different components of the velocity by the Lorentz factor 
and also the coupling between inertial and thermodynamics. We 
demonstrated that Adaptive Mesh Refinement (AMR) is then very 
useful to resolve the associated very thin structures properly. We 
tested the code ability with stringent recent test problems collected 
from the astrophysical literature, including ID and 2D shock tube 
problems, an ultrarelativistic flow reflecting of a wall, a relativistic 
variant of a forward-faced reflecting step, and 2D astrophysical jet 
propagation. 

We used the fireball model to investigate ID and 2D after- 
glow phases in GRBs. In ID, we examined the evolution of a cold 
relativistic shell with a Lorentz factor of 100 in uniform medium. 
In this simulation we discussed details of the internal structure of 
the evolving shell-ISM ejecta and compare them with analytical 
estimates. We followed the evolution of this isotropic explosion al- 
most all the way into the classical Sedov phase. At all times, we 
resolve the various regions that characterize this interaction. We 
quantified and discussed the precise deceleration of the relativistic 
shell. When most of the initial energy of the shell is transferred to 
swept-up shocked ISM (this occurs at the forward shock), the de- 
celeration of shocked ISM is eventually well described by the rela- 
tivistic Blandford-McKee self-similar solution. Hence the Lorentz 
factor of the forward shock decreases as R~ 3/2 . 

We investigate also the afterglow phase for a beamed 2D shell. 
In this model, we discussed analogies and important differences 
with the ID model. The interaction of a confined relativistic shell 
with the ISM is characterised by the appearance of a bow shock. We 
showed how ISM material is laterally pushed out, thus decreasing 
the amount of accumulated matter in front of the shell near the axis. 
The part of the shell furthest away from the axis decelerates then 
faster than in a ID spherical case. Although the deceleration of the 
shell starts early as compared to an equivalent isotropic case, the 
deceleration of the inner part of the shell is slow due to the weak 
accreted ISM matter in front of the shell. The thermal energy of the 
shocked ISM increases and induces a lateral spread of this shocked 
ISM. We have shown with a high resolution simulation of jet-like 
GRB models in their afterglow phase that this lateral expansion 
goes through various phases. 

First, the shell spreads only with its initial lateral velocity un- 
til it accretes enough ISM matter. In this phase, the shocked ISM 
spreads laterally with a velocity near the sound speed. However, 
the reverse shock propagates in a Newtonian fashion through the 
shell, thus having a small efficiency in the conversion of the kinetic 
energy of the shell to thermal energy, hence the expansion of the 
shocked shell is still weak. Only in a later phase when the reverse 
shock becomes relativistic, the lateral expansion of the shocked 
shell increases drastically and reaches a high velocity v T ~ 0.7c. 
The transition from slow to fast lateral spreading of the shell is 
thus related to the transition from Newtonian to relativistic reverse 
shock propagation. However, as the forward shock is always rela- 
tivistic the shocked ISM spreads laterally faster. 

The 2D simulation has revealed rapid internal motions in the 
decelerating configuration. It is possible that these motions result in 
modulations in the afterglow emission. In future work, we intend to 
use these and similar simulation results to compute their predictions 
for the precise afterglow spectral evolution. 
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